home *** CD-ROM | disk | FTP | other *** search
/ Nebula 1 / Nebula One.iso / Graphics / Misc / NeXTcontour_1.7 / Source / contour_plot.f < prev    next >
Encoding:
Text File  |  1995-06-12  |  8.2 KB  |  317 lines
  1.       
  2.       subroutine contour(jdim,kdim,nj,nk,x,y,f,
  3.      $     ncont,acont)
  4.  
  5.       dimension x(jdim,kdim),y(jdim,kdim),f(jdim,kdim)
  6.       dimension acont(ncont)
  7.  
  8.       call con2l(jdim,kdim,1,nj,1,nk,x,y,f,ncont,acont)
  9. c
  10.       return
  11.       end
  12. c
  13.       subroutine cline2(cont,np,x,y)
  14.  
  15.       dimension x(np),y(np)
  16. c
  17. cccc      call concol(cont,funmin,funmax,indexlo,indexup,icolor)
  18. ccc      call move2(x(1),y(1))
  19. ccc          do i = 1,np
  20. ccc           write(3,*)x(i),y(i)
  21. ccc        enddo
  22. ccc 1       call draw2(x(i),y(i))
  23.  
  24.  
  25. c
  26.       return
  27.       end
  28. C
  29. C-----------------------------------------------------------------------
  30.       SUBROUTINE CON2L(IDIM,JDIM,IS,IE,JS,JE,X,Y,F,NCONT,ACONT)
  31. C
  32. C   Calculate contour lines for the function F in the region IS to IE, JS to
  33. C   JE.  X,Y coordinates corresponding to the grid points are in arrays X and
  34. C   Y.
  35. C
  36. C   Common block TEMP1 contains scratch arrays XT, YT, ZT, and IA, and may be
  37. C   used elsewhere.  This subroutine is supposed to be able to recover from
  38. C   filling either array.
  39. C
  40.       PARAMETER (MTEMP1=1000,MIA=3000)
  41.       COMMON /TEMP1/ XT(MTEMP1),YT(MTEMP1),ZT(MTEMP1),IA(MIA)
  42.       DIMENSION X(IDIM,JDIM),Y(IDIM,JDIM),F(IDIM,JDIM)
  43.       DIMENSION ACONT(NCONT)
  44. C
  45. C   One little check.  If IS=IE or JS=JE, return with no contour lines.
  46. C
  47.       IF (IS.EQ.IE .OR. JS.EQ.JE) RETURN
  48. C
  49. C   Zero the marker array.
  50. C
  51.       DO 100 IIA= 1,MIA
  52.          IA(IIA)= 0
  53.   100    CONTINUE
  54. C
  55. C   LNSTRT=1 means we're starting a new line.
  56. C
  57.       LNSTRT= 1
  58.       IGNEXT= 0
  59. C
  60. C   Loop through each contour level.
  61. C
  62.       DO 1000 ICONT= 1,NCONT
  63.          CONT= ACONT(ICONT)
  64. C
  65. C   The IS-IE,JS-JE region may have to be subdivided because of IA space.  This
  66. C   may have a detrimental effect on labelling of the contour lines.  Move 
  67. C   IS,IE,JS,JE to new variables.
  68. C
  69.          ISS= IS
  70.          IEE= IE
  71.          JSS= JS
  72.          JEE= JE
  73. C
  74. C   Flag points in IA where the the function increases through the contour
  75. C   level, not including the boundaries.  This is so we have a list of at least
  76. C   one point on each contour line that doesn't intersect a boundary.
  77. C
  78.   110    CONTINUE
  79.          IAE= 0
  80.          DO 200 J= JSS+1,JEE-1
  81.             IMB  = 0
  82.             IAEND= IAE
  83.             DO 200 I= ISS,IEE
  84.                IF (F(I,J).LE.CONT) THEN
  85.                   IMB    = 1
  86.                ELSE IF (IMB.EQ.1) THEN
  87.                   IAE    = IAE+1
  88.                   IA(IAE)= 1000*I+J
  89.                   IMB    = 0
  90. C
  91. C   Check if the IA array is full.  If so, the subdividing algorithm goes like
  92. C   this:  if we've marked at least one J row, drop back to the last completed
  93. C   J and call that the region.  If we haven't even finished one J row, our 
  94. C   region just extends to this I location.
  95. C
  96.                   IF (IAE.EQ.MIA) THEN
  97.                      IF (J.GT.JSS+1) THEN
  98.                         IAE= IAEND
  99.                         JEE= J
  100.                      ELSE
  101.                         JEE= MIN0(J+1,JEE)
  102.                         IEE= I
  103.                      ENDIF
  104.                      GOTO 210
  105.                   ENDIF
  106. C
  107.                ENDIF
  108.   200          CONTINUE 
  109. C
  110. C   Search along the boundaries for contour line starts.  IMA tells which 
  111. C   boundary of the region we're on.
  112. C
  113.   210    CONTINUE
  114.          IMA = 1
  115.          IMB = 0
  116.          IBEG= ISS-1
  117.          JBEG= JSS
  118. C
  119.     1    IBEG= IBEG+1
  120.          IF (IBEG.EQ.IEE)  IMA= 2 
  121.          GOTO 5
  122.     2    JBEG= JBEG+1
  123.          IF (JBEG.EQ.JEE)  IMA= 3 
  124.          GOTO 5
  125.     3    IBEG= IBEG-1
  126.          IF (IBEG.EQ.ISS)  IMA= 4
  127.          GOTO 5
  128.     4    JBEG= JBEG-1
  129.          IF (JBEG.EQ.JSS)  IMA= 5
  130.     5    IF (F(IBEG,JBEG).GT.CONT) GOTO 7
  131.          IMB= 1
  132.     6    GOTO (1,2,3,4,91),IMA 
  133.     7    IF (IMB.NE.1) GOTO 6
  134. C
  135. C   Got a start point.
  136. C
  137.          IMB= 0
  138.          I  = IBEG 
  139.          J  = JBEG 
  140.          FIJ= F(IBEG,JBEG) 
  141. C
  142. C   Round the corner if necessary.
  143. C
  144.          GOTO (21,11,12,13,51),IMA 
  145.    11    IF (J.NE.JSS) GOTO 31
  146.          GOTO 21 
  147.    12    IF (I.NE.IEE) GOTO 41 
  148.          GOTO 31 
  149.    13    IF (J.NE.JEE) GOTO 51 
  150.          GOTO 41 
  151. C
  152. C   This is how we start in the middle of the region, using IA.
  153. C
  154.    10    I      = IA(IIA)/1000
  155.          J      = IA(IIA)-1000*I  
  156.          FIJ    = F(I,J) 
  157.          IA(IIA)= 0
  158.          GOTO 21 
  159. C                                                                       4
  160. C   Look different directions to see which way the contour line went: 1-|-3
  161. C                                                                       2
  162.    20    J   = J+1
  163.    21    I   = I-1
  164.          IF (I.LT.ISS) GOTO 90
  165.          IDIR= 1
  166.          IF (F(I,J).LE.CONT) GOTO 52 
  167.          FIJ = F(I,J) 
  168.          GOTO 31 
  169.    30    I   = I-1
  170.    31    J   = J-1
  171.          IF (J.LT.JSS) GOTO 90
  172.          IDIR= 2
  173.          IF (F(I,J).LE.CONT) GOTO 60 
  174.          FIJ = F(I,J) 
  175.          GOTO 41 
  176.    40    J   = J-1
  177.    41    I   = I+1
  178.          IF (I.GT.IEE) GOTO 90 
  179.          IDIR= 3
  180.          IF (F(I,J).LE.CONT) GOTO 60 
  181.          FIJ = F(I,J) 
  182.          GOTO 51 
  183.    50    I   = I+1
  184.    51    J   = J+1
  185.          IDIR= 4
  186.          IF (J.GT.JEE) GOTO 90 
  187.          IF (F(I,J).LE.CONT) GOTO 60 
  188.          FIJ = F(I,J) 
  189.          GOTO 21 
  190. C
  191. C   Wipe this point out of IA if it's in the list.
  192. C
  193.    52    IF (IAE.EQ.0) GOTO 60 
  194.          IJ= 1000*I+J+1000 
  195.          DO 300 IIA= 1,IAE 
  196.             IF (IA(IIA).EQ.IJ) THEN
  197.                IA(IIA)= 0
  198.                GOTO 60
  199.             ENDIF
  200.   300       CONTINUE
  201. C
  202. C   Do interpolation for X,Y coordinates.
  203. C
  204.    60    XYF= (CONT-F(I,J))/(FIJ-F(I,J))
  205. C
  206. C   This tests for a contour point coinciding with a grid point.  In this case
  207. C   the contour routine comes up with the same physical coordinate twice.  If
  208. C   If we don't trap it, it can (in some cases significantly) increase the 
  209. C   number of points in a contour line.  Also, if this happens on the first 
  210. C   point in a line, the second point could be misinterpreted as the end of a
  211. C   (circling) contour line.
  212. C
  213.          IF (XYF.EQ.0.) IGNEXT= IGNEXT+1
  214.          GOTO (61,62,63,64),IDIR
  215.    61    WXX= X(I,J)+XYF*(X(I+1,J)-X(I,J))
  216.          WYY= Y(I,J)+XYF*(Y(I+1,J)-Y(I,J))
  217.          GOTO 65 
  218.    62    WXX= X(I,J)+XYF*(X(I,J+1)-X(I,J))
  219.          WYY= Y(I,J)+XYF*(Y(I,J+1)-Y(I,J))
  220.          GOTO 65 
  221.    63    WXX= X(I,J)+XYF*(X(I-1,J)-X(I,J))
  222.          WYY= Y(I,J)+XYF*(Y(I-1,J)-Y(I,J))
  223.          GOTO 65 
  224.    64    WXX= X(I,J)+XYF*(X(I,J-1)-X(I,J))
  225.          WYY= Y(I,J)+XYF*(Y(I,J-1)-Y(I,J))
  226. C
  227. C   Figure out what to do with this point.
  228. C
  229.    65    CONTINUE 
  230.          IF (LNSTRT.NE.1) GOTO 66
  231. C
  232. C   This is the first point in a contour line.
  233. C
  234.          NP    = 1
  235.          XT(NP)= WXX
  236.          YT(NP)= WYY
  237.          WX    = WXX 
  238.          WY    = WYY 
  239.          LNSTRT= 0 
  240.          GOTO 67
  241. C
  242. C   Add a point to this line.  Check for duplicate point first.
  243. C
  244.    66    CONTINUE
  245.          IF (IGNEXT.EQ.2) THEN
  246.             IF (WXX.EQ.XT(NP) .AND. WYY.EQ.YT(NP)) THEN
  247.                IGNEXT= 0
  248.                GOTO 67
  249.             ELSE
  250.                IGNEXT= 1
  251.             ENDIF
  252.          ENDIF
  253. C
  254.          NP    = NP+1
  255.          XT(NP)= WXX
  256.          YT(NP)= WYY
  257. C
  258. C   See if the temporary array xT is full.
  259. C
  260.          IF (NP.EQ.MTEMP1) THEN
  261.             CALL CLINE2(CONT,NP,XT,YT)
  262.             XT(1)= XT(NP)
  263.             YT(1)= YT(NP)
  264.             NP   = 1
  265.          ENDIF
  266. C
  267. C   Check to see if we're back to the initial point.
  268. C
  269.          IF (WXX.NE.WX) GOTO 67
  270.          IF (WYY.EQ.WY) GOTO 90
  271. C
  272. C   Nope.  Search for the next point on this line.
  273. C
  274.    67    CONTINUE
  275.          GOTO (50,20,30,40),IDIR
  276. C
  277. C   This is the end of a contour line.  After this we'll start a new line.
  278. C
  279.    90    LNSTRT= 1
  280.          IGNEXT= 0
  281.          CALL CLINE2(CONT,NP,XT,YT)
  282. C
  283. C   If we're not done looking along the boundaries, go look there some more.
  284. C
  285.          IF (IMA.NE.5) GOTO 6
  286. C
  287. C   Otherwise, get the next start out of IA.
  288. C
  289.    91    CONTINUE
  290.          IF (IAE.NE.0) THEN
  291.             DO 400 IIA= 1,IAE
  292.                IF (IA(IIA).NE.0) GOTO 10
  293.   400          CONTINUE
  294.          ENDIF
  295. C
  296. C   And if there are no more of these, we're done with this region.  If we've
  297. C   subdivided, update the region pointers and go back for more.
  298. C
  299.          IF (IEE.EQ.IE) THEN
  300.             IF (JEE.NE.JE) THEN
  301.                JSS= JEE
  302.                JEE= JE
  303.                GOTO 110
  304.             ENDIF
  305.          ELSE
  306.             ISS= IEE
  307.             IEE= IE
  308.             GOTO 110
  309.          ENDIF
  310. C
  311. C   Loop back for the next contour level.
  312. C
  313.  1000    CONTINUE
  314. C
  315.       RETURN 
  316.       END
  317.